FUT2_Baseline_tissue

Structure and composition

Reference sites

Packages

pacman::p_load(
  dplyr,
  tidyverse,
  readxl,
  vegan)
# setwd("D:/PhD/01_Data/01_Baseline/Model_1") # Windows
# setwd("/Volumes/Yiming_Wang/PhD/01_Data/01_Baseline/Model_1/Tissue") # Mac
# setwd("/Volumes/Yiming_Wang/PhD/01_Data/07_Tissue/Tissue sequencing_control group/All tissue_filtered table") #Imac pro
# setwd("/Users/godric/Desktop/Fut2 review comments/Revision") #Macbook pro
setwd("/Volumes/Yiming_Wang/PhD/15_Manuscript/Submission/07_ISME/Review comments/Revision data/26 April 2024")#Yiming's imac pro
df_all <- read_excel("Baseline_S47_Tissue.xlsx")

df_ST1T <- df_all %>% 
  filter(SampleType == "ST1_T")

df_ST2T <- df_all %>% 
  filter(SampleType == "ST2_T")

df_LT <- df_all %>% 
  filter(SampleType == "LT")

# Subgroup by variables
#30436d WT Male
#9ea5c2 WT Female
#654e3a KO Male
#bcaba6 KO Female

#3C5488B2 WT
#7E6148B2 KO
#D98594FF Female
#94C5CCFF Male

Beta diversity_results and visualisation

Datashape

# places that are not satisfied
  ## 1) Permanova results table
  ## 2) Figure A and B did not align because of legend (Change legend size, position to bottom, change P value position)


#library packages
pacman::p_load(
  data.table,
  vegan,
  knitr,
  EcolUtils
)

# Shape data
# Level orders
# df_all$Genotype <- factor (df_all$Genotype, levels = c("WT", "KO"))
# df_all$Sex <- factor(df_all$Sex, levels = c("Male", "Female"))

df_ST1T$GenotypeSex <- factor(df_ST1T$GenotypeSex, levels = c("WT Male","KO Male"))
df_ST2T$GenotypeSex <- factor(df_ST2T$GenotypeSex, levels = c("WT Male","KO Male"))
df_LT$GenotypeSex <- factor(df_LT$GenotypeSex, levels = c("WT Male","KO Male"))

Bray-curtis

Proximal small intestine

set.seed(123)
# subset dateframe
df_abundance <- subset(df_ST1T, select = -c(1:6))
df_metadata <- subset(df_ST1T, select = c(1:6))

# Create bray-curtis distance matrix
dist_bray <- vegdist(df_abundance,  method = "bray")

# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)

## Run 0 stress 0.04059044 
## Run 1 stress 0.04059042 
## ... New best solution
## ... Procrustes: rmse 8.989373e-05  max resid 0.0001311121 
## ... Similar to previous best
## Run 2 stress 0.0405904 
## ... New best solution
## ... Procrustes: rmse 1.546332e-05  max resid 3.034007e-05 
## ... Similar to previous best
## Run 3 stress 0.2274369 
## Run 4 stress 0.0405904 
## ... Procrustes: rmse 1.309338e-05  max resid 3.302249e-05 
## ... Similar to previous best
## Run 5 stress 0.04059044 
## ... Procrustes: rmse 8.680016e-05  max resid 0.0001219727 
## ... Similar to previous best
## Run 6 stress 0.1139412 
## Run 7 stress 0.3314815 
## Run 8 stress 0.117235 
## Run 9 stress 0.04059063 
## ... Procrustes: rmse 0.0001848809  max resid 0.0002566614 
## ... Similar to previous best
## Run 10 stress 0.117235 
## Run 11 stress 0.04059041 
## ... Procrustes: rmse 2.967561e-05  max resid 4.274677e-05 
## ... Similar to previous best
## Run 12 stress 0.1594032 
## Run 13 stress 0.1139412 
## Run 14 stress 0.04059041 
## ... Procrustes: rmse 3.577167e-05  max resid 7.251582e-05 
## ... Similar to previous best
## Run 15 stress 0.1139412 
## Run 16 stress 0.04059054 
## ... Procrustes: rmse 0.0001563901  max resid 0.0002213148 
## ... Similar to previous best
## Run 17 stress 0.117235 
## Run 18 stress 0.175473 
## Run 19 stress 0.04059043 
## ... Procrustes: rmse 5.695116e-05  max resid 8.059366e-05 
## ... Similar to previous best
## Run 20 stress 0.04059044 
## ... Procrustes: rmse 6.588413e-05  max resid 9.0608e-05 
## ... Similar to previous best
## *** Solution reached

NMDS$stress

## [1] 0.0405904

# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray") 

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
##             Df SumOfSqs      R2      F Pr(>F)
## GenotypeSex  1  0.16453 0.19065 1.8845 0.1252
## Residual     8  0.69844 0.80935              
## Total        9  0.86297 1.00000

# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.006147 0.0061474 0.3078    999  0.504
## Residuals  8 0.159771 0.0199713

Proximal small intestine figure

# Figures
#library packages
pacman::p_load(
  ggthemes,
  ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )

# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
        separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)


centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
      separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)

seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
             by = c('GenotypeSex','Genotype','Sex'), sort = FALSE) 


composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
  geom_segment(data = seg,
             mapping = aes(xend = oNMDS1, yend = oNMDS2),
             size = 0.25) +
  geom_point(data = centroids, size = 4) +
  geom_point()+ #geo,_point()
  stat_ellipse(linetype="longdash",size=0.3,color="gray32")+
  coord_fixed()+
  theme_bw() +
  labs(title="")+
  theme(legend.position = "bottom",
        legend.title= element_text(size=10, face="plain"),
        legend.text = element_text(size=14),
        legend.key = element_rect(fill = "transparent", colour = "black"),
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 9), size=16),
        axis.title.y=element_text(margin = margin(r = 3), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_x_continuous(limits = c(-1.5,1.5),breaks = seq(-1.5,1.5,0.6))+
  scale_y_continuous(limits = c(-2.0,2.0),breaks = seq(-2.0,2.0,0.8))+
  scale_color_manual(values=c("#30436d","#654e3a"))+
  scale_fill_manual(values=c("#30436d","#654e3a"))+
  scale_shape_manual(values=c(24,21,24,21))+
  annotate("text",x=-1.5,y=-2,
           label="Stress = 0.041",hjust=0, size=4)+
  annotate("text",x=-1.5,y=2,
           label="Proximal small intestine",hjust=0, size=5, fontface = "bold")+
  annotate("text",x=-1.5,y=1.6,
           label=expression(paste("Male (WT vs KO): ",italic("P") > 0.05)),hjust=0, size=4)

# Print
composition_GenotypeSex_final

Distal small intestine

set.seed(123)
# subset dateframe
df_abundance <- subset(df_ST2T, select = -c(1:6))
df_metadata <- subset(df_ST2T, select = c(1:6))

# Create bray-curtis distance matrix
dist_bray <- vegdist(df_abundance,  method = "bray")

# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)

## Run 0 stress 0.04040786 
## Run 1 stress 0.04040786 
## ... New best solution
## ... Procrustes: rmse 8.780053e-06  max resid 1.891585e-05 
## ... Similar to previous best
## Run 2 stress 0.1231312 
## Run 3 stress 0.04040786 
## ... New best solution
## ... Procrustes: rmse 3.686835e-06  max resid 7.196533e-06 
## ... Similar to previous best
## Run 4 stress 0.04040786 
## ... Procrustes: rmse 7.393811e-06  max resid 1.550481e-05 
## ... Similar to previous best
## Run 5 stress 0.04040786 
## ... Procrustes: rmse 8.841731e-06  max resid 1.950637e-05 
## ... Similar to previous best
## Run 6 stress 0.04040786 
## ... New best solution
## ... Procrustes: rmse 2.12053e-06  max resid 4.289854e-06 
## ... Similar to previous best
## Run 7 stress 0.1231312 
## Run 8 stress 0.04040787 
## ... Procrustes: rmse 3.053269e-05  max resid 6.377193e-05 
## ... Similar to previous best
## Run 9 stress 0.0472167 
## Run 10 stress 0.1252416 
## Run 11 stress 0.04040786 
## ... Procrustes: rmse 1.870057e-06  max resid 4.302461e-06 
## ... Similar to previous best
## Run 12 stress 0.04040786 
## ... Procrustes: rmse 2.486584e-06  max resid 4.60433e-06 
## ... Similar to previous best
## Run 13 stress 0.04040786 
## ... Procrustes: rmse 2.887797e-06  max resid 5.481577e-06 
## ... Similar to previous best
## Run 14 stress 0.04040786 
## ... Procrustes: rmse 9.520293e-06  max resid 1.977447e-05 
## ... Similar to previous best
## Run 15 stress 0.04040786 
## ... Procrustes: rmse 2.542568e-06  max resid 4.669764e-06 
## ... Similar to previous best
## Run 16 stress 0.0472167 
## Run 17 stress 0.04040786 
## ... Procrustes: rmse 3.249232e-06  max resid 5.789046e-06 
## ... Similar to previous best
## Run 18 stress 0.1231313 
## Run 19 stress 0.04040787 
## ... Procrustes: rmse 1.450243e-05  max resid 2.672569e-05 
## ... Similar to previous best
## Run 20 stress 0.1231313 
## *** Solution reached

NMDS$stress

## [1] 0.04040786

# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray") 

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
##             Df SumOfSqs      R2     F Pr(>F)   
## GenotypeSex  1  0.58745 0.54257 9.489 0.0071 **
## Residual     8  0.49527 0.45743                
## Total        9  1.08272 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
## Groups     1 0.030943 0.0309426 4.2672    999  0.082 .
## Residuals  8 0.058010 0.0072512                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Distal small intestine figure

# Figures
#library packages
pacman::p_load(
  ggthemes,
  ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )

# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
        separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)


centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
      separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)

seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
             by = c('GenotypeSex','Genotype','Sex'), sort = FALSE) 


composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
  geom_segment(data = seg,
             mapping = aes(xend = oNMDS1, yend = oNMDS2),
             size = 0.25) +
  geom_point(data = centroids, size = 4) +
  geom_point()+ #geo,_point()
  stat_ellipse(linetype="longdash",size=0.3,color="gray32")+
  coord_fixed()+
  theme_bw() +
  labs(title="")+
  theme(legend.position = "bottom",
        legend.title= element_text(size=10, face="plain"),
        legend.text = element_text(size=14),
        legend.key = element_rect(fill = "transparent", colour = "black"),
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 9), size=16),
        axis.title.y=element_text(margin = margin(r = 3), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_x_continuous(limits = c(-1.5,1.5),breaks = seq(-1.5,1.5,0.6))+
  scale_y_continuous(limits = c(-2,2),breaks = seq(-2,2,0.8))+
  scale_color_manual(values=c("#30436d","#654e3a"))+
  scale_fill_manual(values=c("#30436d","#654e3a"))+
  scale_shape_manual(values=c(24,21,24,21))+
  annotate("text",x=-1.5,y=-2,
           label="Stress = 0.040",hjust=0, size=4)+
  annotate("text",x=-1.5,y=2,
           label="Distal small intestine",,hjust=0, size=5, fontface = "bold")+
  annotate("text",x=-1.5,y=1.6,
           label=expression(paste("Male (WT vs KO): ",italic("P") == 0.0071)),hjust=0, size=4)

# Print
composition_GenotypeSex_final

Large intestine

set.seed(123)
# subset dateframe
df_abundance <- subset(df_LT, select = -c(1:6))
df_metadata <- subset(df_LT, select = c(1:6))
 
# Create bray-curtis distance matrix
dist_bray <- vegdist(df_abundance,  method = "bray")


# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)

## Run 0 stress 0.119291 
## Run 1 stress 0.119291 
## ... Procrustes: rmse 9.908197e-05  max resid 0.0001913734 
## ... Similar to previous best
## Run 2 stress 0.119291 
## ... New best solution
## ... Procrustes: rmse 3.670601e-05  max resid 5.712139e-05 
## ... Similar to previous best
## Run 3 stress 0.2385138 
## Run 4 stress 0.119291 
## ... Procrustes: rmse 0.0001347919  max resid 0.0002617218 
## ... Similar to previous best
## Run 5 stress 0.119291 
## ... Procrustes: rmse 1.54859e-05  max resid 2.286731e-05 
## ... Similar to previous best
## Run 6 stress 0.2418937 
## Run 7 stress 0.1192911 
## ... Procrustes: rmse 6.660418e-05  max resid 0.0001270189 
## ... Similar to previous best
## Run 8 stress 0.2825862 
## Run 9 stress 0.119291 
## ... Procrustes: rmse 4.138662e-07  max resid 7.575965e-07 
## ... Similar to previous best
## Run 10 stress 0.1962946 
## Run 11 stress 0.1192915 
## ... Procrustes: rmse 0.0001413468  max resid 0.0002531473 
## ... Similar to previous best
## Run 12 stress 0.2444607 
## Run 13 stress 0.2532164 
## Run 14 stress 0.119291 
## ... Procrustes: rmse 6.078828e-05  max resid 0.0001113746 
## ... Similar to previous best
## Run 15 stress 0.119291 
## ... Procrustes: rmse 7.885862e-05  max resid 0.0001437482 
## ... Similar to previous best
## Run 16 stress 0.119291 
## ... Procrustes: rmse 7.812125e-05  max resid 0.0001434385 
## ... Similar to previous best
## Run 17 stress 0.119291 
## ... Procrustes: rmse 6.768069e-06  max resid 1.245338e-05 
## ... Similar to previous best
## Run 18 stress 0.119291 
## ... Procrustes: rmse 3.591897e-05  max resid 5.235938e-05 
## ... Similar to previous best
## Run 19 stress 0.2518441 
## Run 20 stress 0.119291 
## ... New best solution
## ... Procrustes: rmse 1.59648e-05  max resid 3.229225e-05 
## ... Similar to previous best
## *** Solution reached

NMDS$stress

## [1] 0.119291

# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray") 

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
##             Df SumOfSqs     R2      F Pr(>F)
## GenotypeSex  1  0.12799 0.1625 1.5522 0.1793
## Residual     8  0.65962 0.8375              
## Total        9  0.78761 1.0000

# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)

## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.0000402 0.00004020 0.0421    999  0.921
## Residuals  8 0.0076406 0.00095507

Large intestine figure

# Figures
#library packages
pacman::p_load(
  ggthemes,
  ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )

# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
        separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)


centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
      separate(GenotypeSex,c("Genotype","Sex"),sep=" ",remove=FALSE)

seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
             by = c('GenotypeSex','Genotype','Sex'), sort = FALSE) 


composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
  geom_segment(data = seg,
             mapping = aes(xend = oNMDS1, yend = oNMDS2),
             size = 0.25) +
  geom_point(data = centroids, size = 4) +
  geom_point()+ #geo,_point()
  stat_ellipse(linetype="longdash",size=0.3,color="gray32")+
  coord_fixed()+
  theme_bw() +
  labs(title="")+
  theme(legend.position = "bottom",
        legend.title= element_text(size=10, face="plain"),
        legend.text = element_text(size=14),
        legend.key = element_rect(fill = "transparent", colour = "black"),
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 9), size=16),
        axis.title.y=element_text(margin = margin(r = 3), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_x_continuous(limits = c(-1.5,1.5),breaks = seq(-1.5,1.5,0.6))+
  scale_y_continuous(limits = c(-2,2),breaks = seq(-2,2,0.8))+
  scale_color_manual(values=c("#30436d","#654e3a"))+
  scale_fill_manual(values=c("#30436d","#654e3a"))+
  scale_shape_manual(values=c(24,21,24,21))+
  annotate("text",x=-1.5,y=-2,
           label="Stress = 0.12",hjust=0, size=4)+
  annotate("text",x=-1.5,y=2,
           label="Large intestine",,hjust=0, size=5, fontface = "bold")+
  annotate("text",x=-1.5,y=1.6,
           label=expression(paste("Male (WT vs KO): ",italic("P") > 0.05)),hjust=0, size=4)

# Print
composition_GenotypeSex_final

ANOISM: Variance between groups and within groups

• ANOSIM analysis to test whether the variance between groups is larger than within groups (if it is larger, then P values is less than 0.05)
• Anoism analysis (R value from -1 to 1, R > 0 indicates between-group difference > within group difference; R < 0 indicates within group > between group difference; P value indicates the statistical significance)

# ST1T
## Define group
df_abundance <- subset(df_ST1T, select = -c(1:6))
df_metadata <- subset(df_ST1T, select = c(1:6))

## Bray-curtis
dist_bray <- vegdist(df_abundance,  method = "bray")

## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)

## 
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R:   0.1 
##       Significance: 0.174 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.184 0.252 0.324 0.332 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%  75% 100%  N
## Between  1 14.00 26.0 34.0   45 25
## WT Male  3 18.25 22.0 34.5   44 10
## KO Male  7  9.25 14.5 29.5   37 10

plot(anosim.GenotypeSex_bray)

# ST2T
## Define group
df_abundance <- subset(df_ST2T, select = -c(1:6))
df_metadata <- subset(df_ST2T, select = c(1:6))

## Bray-curtis
dist_bray <- vegdist(df_abundance,  method = "bray")

## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)

## 
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.812 
##       Significance: 0.007 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.236 0.288 0.408 0.460 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%  75% 100%  N
## Between 16 25.00 31.0 39.0   45 25
## WT Male  1  4.25  6.5  9.5   14 10
## KO Male  3 11.50 18.0 23.5   37 10

plot(anosim.GenotypeSex_bray)

# LT
## Define group
df_abundance <- subset(df_LT, select = -c(1:6))
df_metadata <- subset(df_LT, select = c(1:6))

## Bray-curtis
dist_bray <- vegdist(df_abundance,  method = "bray")

## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)

## 
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.12 
##       Significance: 0.768 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.203 0.284 0.308 0.392 
## 
## Dissimilarity ranks between and within classes:
##         0%   25% 50%   75% 100%  N
## Between  1  9.00  18 39.00   45 25
## WT Male  3 17.75  29 34.75   38 10
## KO Male  8 21.25  24 28.50   33 10

plot(anosim.GenotypeSex_bray)